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The last two decades have witnessed tremendous growth in computational power, the 
development of computational fluid dynamics (CFD) codes which scale well over thousands 
of processors, and the refinement of unstructured grid-generation tools which facilitate 
rapid surface and volume gridding of complex geometries. Thus, engineering calculations 
of 10 7 — 10 8 finite-volume cells have become routine for some types of problems. Although 
the Reynolds Averaged Navier Stokes (RANS) approach to modeling turbulence is still in 
extensive and wide use, increasingly large-eddy simulation (LES) and hybrid RANS-LES 
approaches are being applied to resolve the largest scales of turbulence in many engineering 
problems. However, it has also become evident that LES places different requirements on 
the numerical approaches for both the spatial and temporal discretization of the Navier 
Stokes equations than does RANS. In particular, LES requires high time accuracy and 
minimal intrinsic numerical dispersion and dissipation over a wide spectral range. In this 
paper, the performance of both central-difference and upwind-biased spatial discretizations 
is examined for a one-dimensional acoustic standing wave problem, the Taylor-Green vortex 
problem, and the turbulent channel flow problem. 


Nomenclature 


a sound speed 

E specific total energy, e + 

Ek kinetic energy, U 2 / 2 

e specific internal energy 

f body force per unit volume vector 

h channel half- height, reference length for channel flow problem 
k wavenumber 

l reference length for Taylor- Green vortex problem 

Mach number, U/a 
mass flow rate 

Number of grid points in a particular direction 
p pressure 

Q Q-criterion, (R : R — S : S)/2 
q heat flux 

R gas constant 

R rotation-rate tensor, (Vu — Vu T )/2 

S strain-rate tensor, (Vu + Vu T )/2 

T temperature 

t time 

U velocity magnitude, ||u|| 

u velocity vector 
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V volume 

u,v,w components of velocity vector in the x-, y -, and ^-direction, respectively 
x,y,z cartesian coordinates 

Subscripts 

0 denotes an initial condition 

c denotes a flow property at the centerline 

i, j, k computational indices in the x-, y-, and z-direction, respectively 
m denotes a mean flow property 
w denotes a flow property at the wall 

Symbols 

S kronecker delta 

7 ratio of specific heats 

A thermal conductivity 

fi viscosity 

Q enstrophy, uj • uj/2 

l o vorticity vector 

0 generic flowfield quantity 

p density 

r viscous stress tensor 

5 turbulent kinetic energy dissipation rate 

Superscripts 

* variable normalized by a characteristic flow timescale 

+ variable normalized by inner- law variables 

I. Introduction 

T urbulence is a critical factor in most aero- and propulsion-related flows. For many years, the industry- 
standard approach to incorporating the effect of turbulence in computational fluid dynamics (CFD) 
calculations has been the Reynolds Averaged Navier Stokes (RANS) modeling approach. RANS implicitly 
time- averages turbulent motion, and models its effect on the mixing of species, momentum and energy 
in a flow field. The approach is economical and often allows a CFD calculation to proceed to steady- 
state. However, it makes several simplifying assumptions, and the results will only be as good as the 
RANS turbulence model. Meanwhile, the last two decades have witnessed rapid and explosive growth in 
computational power, the development of CFD codes which scale well over thousands of processors, and 
the refinement of unstructured grid-generation tools which facilitate rapid surface and volume gridding of 
complex geometries. CFD engineering calculations of 10 7 - 10 8 finite- volume cells have become routine. 
Thus, large eddy simulation (LES), which attempts to directly resolve the unsteady motion of the largest 
scales of turbulence, is increasingly of interest for many fluids engineering problems. 1 3 Additionally, hybrid 
RANS-LES approaches such as Spalart’s Detached Eddy Simulation (DES), 4 which attempt to marry the 
best strengths of both methods, are in active use. 

Because it seeks to accurately resolve the unsteady motion of the largest scales of turbulence, LES places 
certain requirements on the numerical method used to discretize and advance the Navier Stokes equations. 
In particular, LES requires high time accuracy and minimal intrinsic numerical dispersion and dissipation 
over a wide range of length scales. A number of studies in the literature 5 10 have investigated the suitability 
of various upwind-based shock-capturing schemes for LES. The purpose of this paper is to evaluate the 
suitability of several different central difference and upwind-biased schemes for a simple one-dimensional 
acoustic standing wave problem, as well as DNS and LES of the Taylor-Green vortex problem and a turbulent 
channel flow problem. 
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II. Governing Equations 


The compressible Navier-Stokes equations for a single chemical species are the governing 
study. The time-dependent differential form of these equations is written as follows: 

model for this 


t +V >u)=° 

(i) 


- + V • (puu + pS) = V • t + f 

(2) 


Q t + v • (pu(E + £>)) = V • ( u • r q) + f • u 

(3) 

A calorically perfect ideal gas equation of state is used: p = pRT, and T = (7 — 1 )e/R, 
R = 287 J/Kg-K for all cases studied in this work. The viscous stress tensor, r is defined as 

. 7 = 1.4 and 


r = 2/iS — -/z(V • u)5 

0 

(4) 

and the heat flux vector as 

q = -AVT 

(5) 


The transport properties p and A are calculated using the Sutherland formulas for air. 


III. Numerical Method 

The Navier-Stokes equations are solved by the finite-difference method in the conventional cartesian 
coordinate system. The inviscid flux (the second term on the left-hand side of Eqs. 1, 2, and 3) is the portion 
of the equations that deals with convection of mass, momentum and energy. Historically, the discretization 
of these terms has been a driving issue in CFD research since its inception. Four central difference (CD) 
schemes, four upwind-biased (UB) schemes, and the Fromm scheme for the inviscid flux are evaluated in 
this work. The stencil coefficients for these schemes are shown in table 1 for a generic flowfield quantity, <f > . 
The number in the name for each scheme denotes the formal order of accuracy. Conservative differencing 
was used such that, at node z , j , /c , the derivative in the ^-direction is given by: 

d(t)jj,k _ <t>i+l/2j,k ~ 4>i-l/2,j,k , x 

dx Ax 

where (t>i+i/2j,k an d 4>i-i/2,j,k are interpolated and formed in such a way that the final difference is 
equivalent to one of the stencils shown in Table 1. For example, </>i+i/ 2 ,y/c = (0i,j,/c + $i+ij,k)/2 and 
4>i-i/2,j,k = (0i-i,j,/c + ^j,/c)/2 yields the CD -2 formula, while 4> i+ i/ 2 j,k = (-(t>i-i,j,k + 5<t>i,j,k + 2</> i+1 j i k)/6 
and 4>i-i/2j,k — (— « 4>i-2,j,k + 5^-qy/c + 2^ J? / c )/6 yields the UB-3 formula. It should be noted that the UB 
and Fromm schemes shown assume convection in the positive- x direction. A reflected stencil (about 4>i : j,k) 
would be employed for convection in the negative- x direction. 

Note that the inviscid flux terms contain combinations of several different primitive flow variables. For 
example, in the x-direction, the terms to be differenced are 0 = [pu, pu 2 + p, puv, puw, pu(E + p)] T - In 
divergence form central differencing, the flux values themselves are directly differenced. In skew-symmetric 
form central differencing, the flux is formed from a combination of central differences of the flux values and the 
component flow property variables themselves. The central difference schemes investigated here employ the 
skew-symmetric approach. This approach has been investigated by several researchers in recent years , 11-16 
and appears to be more accurate and stable than divergence form central differencing. In particular, the 
implementation of Ref. 15 for the skew- symmetric scheme of Ref. 13 was used. Additional detail will be 
provided in the final paper. 

Upwinding can be implemented by various methods in a CFD code . 17, 18 In this work, first a leftward- 
biased interpolation for the primitive flow variables = p, u, v , re, and p at the half- node is formed which 
is consistent with the UB and Fromm formulas from Table 1. The remaining state variables are calculated 
using the perfect gas law, and the inviscid flux terms are then calculated from these components. The process 
is then repeated for a rightward-biased interpolation for a half-node using a reflected stencil. The final value 
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Table 1. Coefficients for central difference and upwind-biased stencils for dcfri j^/dx = (c Pi+i/2,j,k ~ 4>i— 1 / 2 , j,k)/^ x - The 
upwind-biased (UB) and Fromm stencils assume convection in the positive-a? direction. 


Name 

<\>i - 4 

4>i - 3 

0Z-2 

4>i - 1 

(pi 

0i+l 

(pi + 2 

0i+3 

4 

CD-2 




-1/2 


1/2 




CD-4 



1/12 

-2/3 


2/3 

-1/12 



CD-6 


-1/60 

3/20 

-3/4 


3/4 

-3/20 

1/60 


CD-8 

1/280 

-4/105 

1/5 

-4/5 


4/5 

-1/5 

4/105 

-1/280 

UB-1 




-1 

1 





TJB-3 



1/6 

-1 

1/2 

1/3 




UB-5 


-1/30 

1/4 

-1 

1/3 

1/2 

-1/20 



UB-T 

1/140 

-7/105 

3/10 

-1 

1/4 

3/5 

-1/10 

1/105 


Fromm- 2 



1/4 

-5/4 

3/4 

1/4 







(a) Dispersion Error (b) Dissipation Error 

Figure 1. Fourier analysis of spatial discretization schemes studied in this work. 


of the inviscid flux at the half- node is provided by the Roe approximate Riemann solver, 19 using these two 
separate left and right states. Again, additional detail on the method will be provided in the final paper. 

Plots of the dispersion error and dissipation error from a fourier analysis of the various schemes are shown 
in figure 1. The scaled wavenumber is equal to 2i rk/N, such that it is equal to n at the spatial Nyquist 
frequency. As discussed by Li, 20 the dispersion error of each of the UB schemes is equal to that of the CD 
scheme of the next higher order of accuracy. All of the CD schemes are free of dissipation error, while the 
UB and Fromm schemes experience an increasing dissipation error at higher wavenumbers. The higher-order 
UB schemes have a significantly smaller dissipation error than the lower-order ones. 

The viscous terms were discretized using 2 nd order central differences. All schemes were advanced using 
a low-storage 4-stage, 2 nd order accurate Runge-Kutta time advancement algorithm. 
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IV. Results and Discussion 


IV. A. One-Dimensional Acoustic Standing Wave Problem 

As a simple means of evaluating the ability of the various inviscid flux schemes to resolve high-frequency 
acoustics or flowfield variations, they were tested on a one-dimensional acoustic standing wave problem. 
This problem assumes a one-dimensional computational domain with periodic boundary conditions. 128 
grid points are used. The domain is initialized with air at p = 0.25 atm, T = 298.15 K, and a small 
sinusoidal velocity variation as a function of axial distance, x, and specified wavenumber, k: u(x, t = 0) = 
(0.1 m/s) cos (kx). Under the assumption of isentropic, small disturbance inviscid flow, the Navier-Stokes 
equations can be simplified to the linearized equations of gas dynamics, and an exact solution derived: 
u(x,t) = (0.05 m/s)(cos(&(x — at)) -\-cos(k(x-\-at))) = (O.lm/s) cos (kx) cos (/cat), where t is time and a is the 
sound speed (346.117 m/s). Due to space limitations in this abstract, the results from this analysis cannot 
be shown. However, in summary, the results clearly demonstrate the fourier characteristics shown in figure 1. 
The exact solution can even be modified to account for the dispersion and dissipation error of a particular 
scheme, and the resulting equation has excellent agreement with the numerical results. The complete results 
from this analysis will be presented in the final paper. 


IV.B. Taylor-Green Vortex Problem 

The Taylor-Green vortex problem is a benchmark case which, from a smooth initial condition, simulates 
vortex stretching and transition to turbulence. It was used as a test case for the 1 st International Workshop 
on High-Order CFD Methods, held at the 50 th AIAA Aerospace Sciences Meeting, January 7-8, 2012, in 
Nashville, Tennessee. DNS results have been calculated by Brachet et al. 21 and, more recently by van Rees 
et al. 22 The latter results are used as the basis for comparison in this paper. 

The solution is computed on a periodic cube domain, — ttI < x, y, z < nl. The flowfield is initialized with: 


u 

v 

w 

V 


'/““(fMIMf). 

— f/„COS (y) sin (y) COS (y) , 
0 , 


Po + 


PoUl 

16 




( 7 ) 


where here l = 0.01m, p 0 = 7271 Pa, T 0 = 298.15 K and Uo = 34.612 m/s. The flow is effectively incom- 
pressible, Mo = Uo/ao = 0.1. The initial condition for p is calculated by the perfect gas law, assuming fixed 
T = Tq. The transport properties are fixed at their values at Tq throughout the calculation. The resultant 
Reynolds number is Re = poU^l/po = 1600. 

The time evolution of the Taylor- Green vortex flowfield can be visualized by showing iso-surfaces of zero 
Q-criterion ( Q = 0) at various times (figure 2). Q > 0 indicates regions of a flowfield in which vorticity 
dominates over strain. This calculation was performed using the CD-4 scheme on a 192 3 point grid, with 
uniform grid spacing throughout the domain. Time is normalized by a characteristic timescale, t* = tUo/l , 
and the iso-surfaces are colored with the normalized velocity magnitude, U/Uo. It is evident from the figure 
that the initial large-scale (k = 1) vortex structure at t* = 0 steadily evolves toward smaller vortex structures 
(t* =4). By t* = 8, the flow has transitioned to turbulence, and proceeds steadily to a more isotropic state. 
The symmetries in the flowfield, discussed by Ref. 21, are evident from figure 2. Also note that the reduction 
in overall kinetic energy as the simulation progresses is visible in the velocity color scale . 

Time histories for the mean kinetic energy, kinetic energy dissipation rate, and enstrophy for the CD-4 
and UB-3 schemes on a 192 3 point grid are compared with the spectral DNS (on a 512 3 point grid) of Ref. 22 
in figures 3, 4, and 5, respectively. The mean kinetic energy for the domain is given by 


1 f U 2 

Ek ' m = ~poV J P T dV 


(8) 


and the mean kinetic energy dissipation rate is £ = The mean enstrophy for the domain is given 
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(d) t* = 12 (e) t* = 16 (f) t* = 20 

Figure 2. Iso-surfaces of zero Q-criterion showing the time evolution of the Taylor-Green vortex problem at Re = 1600, 
computed with the CD-4 scheme on a 192 3 grid. Iso-surfaces are colored with the normalized velocity magnitude, 
U/Uq. t* is a normalized timescale, tUo/l 



(a) Overview (b) Closeup of t* = 7.5 to t* = 10 

Figure 3. Plots showing the time evolution of the normalized mean kinetic energy, Ek,m/^o for the reference spectral 
DNS on a 512 3 point grid, and the CD-4 and UB-3 schemes on 192 3 point grids. 


by 




i 

ToV 


J^V 


(9) 
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(a) Overview (b) Closeup of t* = 7.5 to t* = 10 

Figure 4. Plots showing the time evolution of the normalized mean kinetic energy dissipation rate, eI/Vq for the 
reference spectral DNS on a 512 3 point grid, and the CD-4 and UB-3 schemes on 192 3 point grids. 



(a) Overview (b) Closeup of t* = 7.5 to t* = 10 

Figure 5. Plots showing the normalized mean enstrophy, Cll 2 /Vq for the reference spectral DNS on a 512 3 point grid, 
and the CD-4 and UB-3 schemes on 192 3 point grids. 


It is evident from figures 3, 4, and 5 that the CD-4 scheme at this grid resolution is in very good agreement 
with the reference DNS results for all three mean flow properties. The UB-3 scheme follows the mean kinetic 
energy history reasonably well, although it consistently dissipates the kinetic energy faster than the DNS 
through t* = 8.6. The mean enstrophy history for the UB-3 scheme begins to diverge below the DNS at 
t* =3.5, and remains significantly smaller throughout the remainder of the simulation. This is likely due to 
the dissipative nature of the UB-3 scheme at higher wavenumbers, which become more and more significant 
as the simulation progresses. 

Thus far, simulations for all schemes have been completed on a 128 3 point grid, the CD-4, CD-6 and 
UB-3 schemes on a 192 3 point grid, and the CD-6 scheme on a 256 3 point grid. In the final paper, the 
schemes will be compared on grids ranging from 64 3 - 256 3 points to assess the suitability of the schemes 
for both DNS and LES of the Taylor- Green vortex. 

IV. C. Turbulent Channel Flow Problem 

Turbulent channel flow has been extensively studied by experiments 23 and DNS studies, 24,25 and the mean 
flow profile and statistics have been very well characterized. Therefore, it is an excellent benchmark case for 
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Figure 6. Iso-surfaces of zero Q-criterion showing turbulent channel flow at Re = 14300, computed with the CD-2 
scheme on a 128 X 129 X 128 point grid, a) initial condition, b) statistically stationary state. Iso-surfaces are colored 
with the normalized velocity magnitude, U/um. t* is a normalized timescale, t*tu rn /27vh, representing the number of 
streamwise flowthrough times. 


evaluating various inviscid flux schemes for DNS or LES of wall-bounded shear flows. Here, turbulent channel 
flow of air at p = 0.25 atm, T = 298.15 K and a mean streamwise velocity u m = 44.44 m/s is simulated. 
The flow is essentially incompressible (M m = 0.127). The computational domain is 0 < x < 2tt ft in the 
streamwise (pc-) direction, — ft < y < ft in the vertical (y-) direction, and 0 < z < irh in the spanwise (z-) 
direction. The channel half-width, ft, is 0.01 m. The Reynolds number based on the full channel width, 2 ft, 
is Re = p rn u rn 2h/ p rn — 14300. 

Periodic boundary conditions are used in both the x- and ^-directions. Adiabatic, no-slip wall boundary 
conditions are imposed on the bottom and top (y~) surfaces. In order to prevent the viscous losses at both 
walls from slowly “wearing down” the flow, the integrated area-sum of these forces is computed at each time 
step, and added in as a body force term, /, over the entire computational volume. This term effectively 
reproduces the effect of the pressure gradient needed to sustain a real turbulent channel flow, and maintains 
a constant mass flow rate in the simulation. 

The flow is initialized using a power law mean velocity profile, along with unit wavenumber (k = 1) 
oscillations in y and z, and k = 2 wavenumber oscillations in x. This divergence-free initial condition was 
borrowed from Moin and Kim. 26 


u 

v 

w 



( 10 ) 


where here the initial centerline streamwise velocity, u c ^ o = 50 m/s and Uo = 5 m/s. 

This initial condition, as well as a fully developed turbulent channel flow, can be visualized by showing 
iso-surfaces of zero Q-criterion (Q = 0) at t* = 0 and t* = 30 (figure 6). t* is a normalized timescale, 
t* = ftu m /27rft, representing the number of streamwise flowthrough times. The iso-surfaces are colored with 
the normalized velocity magnitude, U/u m . This calculation was performed using the CD-2 scheme on a 
128 x 129 x 128 point grid. The grid spacing is uniform in x and z, and stretched by the hyperbolic tangent 
function in y. The wall- normal grid spacing at the walls is 0.001ft. Grid spacings normalized by the wall 
friction lengthscale are: Ax + = A Xy/lr w p w ) / fi w = 20.1, A y+ = 0.426, A y+ = 16.4, Az + = 10.0. 

Skin friction coefficient histories for the various inviscid flux schemes are shown in figure 7. The CD 
schemes were started from the initial condition, and it is evident from the skin friction history that the flow 
rapidly transitions to turbulence within the first few flowthrough times. The solutions were time-averaged 
from t* = 15-30, as well as space- averaged in the periodic directions x and z. The time- and space- averaged 
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Figure 7. Comparison of skin friction coefficient history for 




the inviscid flux schemes on a 128 X 129 X 128 point grid. 



Figure 8. Comparison of streamwise velocity profiles for the inviscid flux schemes on a 128 X 129 X 128 point grid. 


skin friction coefficient for all three CD schemes is in excellent agreement (within 3%) with the predictions 
from the reference Dean 23 correlation at this Reynolds number. The UB and Fromm schemes were started 
from the CD-2 solution at t* = 10. As might be expected, the UB-1 scheme rapidly begins to laminarize 
the flow, but has not achieved a steady state condition even by t* = 30. The UB-3, UB-5 and UB-7 
schemes do achieve a stationary state by t* = 15, and their time- and space- averaged coefficient values are 
in progressively closer agreement with the CD-2 results and the prediction from the Dean correlation. The 
average skin friction for the Fromm-2 scheme is slightly lower than the UB-3 scheme. 

The time- and space- averaged streamwise velocity profiles for the schemes (except UB-1 which is only 
space-averaged at t* = 30) are compared with the laminar sublayer profile, the log-law profile, and the DNS 
results of Moser et al. 25 in figure 8. The velocity profiles are normalized by the characteristic wall friction 
velocity, u + = u/y^{T w /p w ), and the distance from the wall is normalized by the wall friction lengthscale, 

= \y ~ yw\\fijwpw) / The results are consistent with the observations from the skin friction history. 

Thus far, simulations for all of schemes have been completed on a 128 x 129 x 128 point grid, the CD-2 
scheme on a 192 x 193 x 192 point grid, and the CD-2, CD-4 and CD-6 schemes on progressively coarser 
96 x 129 x 96, 64 x 129 x 64, and 48 x 129 x 48 point grids. In the final paper, the suitability of all of the 
schemes for DNS and LES of turbulent channel flow will be assessed. The effect of the Smagorinsky and 
Dynamic Smagorinsky subgrid models on the results will also be determined. 
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